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Finding the sources of sound in large nonlinear fields via direct simulation currently requires 
excessive computational cost. This paper describes a simple technique for efficiently solving the 
multidimensional nonlinear Euler equations that significantly reduces this cost and demonstrates 
a useful approach for validating high order nonlinear methods. Up to 15** order accuracy in space 
and time methods were compared and it is shown that an algorithm with a fixed design accuracy 
approaches its maximal utility and then its usefulness exponentially decays unless higher accuracy 
is used. It is concluded that at least a 7 th order method is required to efficiently propagate a 
harmonic wave using the nonlinear Euler equations to a distance of 5 wavelengths while maintaining 
an overall error tolerance that is low enough to capture both the mean flow and the acoustics. 

1. Introduction 

It is well known that aeroacoustic computations require high order numerical schemes to 
efficiently capture the multiscale phenomena occurring over a relatively large time interval. 
But the question of how high has not been answered because of a number of reasons, such 
as the following: 

• Acoustics may be mathematically solved in a number of ways such as via acoustic 
analogies \ direct simulations 2 , and perturbation techniques 3 . 

• Commonly used numerical approaches such as Compact 4 , Essentially Non- Oscillatory 5 , 
Dispersion Relation Preserving 6 , and Galerkin 7 schemes are difficult to extend to 
very high order. 

• The numerical requirements vary considerably due to the wide variety of acoustical 
phenomena which include linear propagation and scattering, inviscid gust interactions 
with objects, receptivity effects, and turbulent noise generation. 

Since most Computational Aeroacou stic s (CAA) techniques have a fixed design accuracy 
and resolution, they offer limited opportunities for exploring the effects of design accuracy on 
the simulation. Instead, an arbitrary accuracy Modified Expansion Solution Approximation 
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(MESA) technique 8,9,10 is used here. This method was recently implemented 11 in a way 
that greatly simplified its use and demonstrated the possibility of subgrid scale resolution 
by including the spatial derivatives of the primitive variables on the grid and by using 
very high design accuracies 12 . For these reasons, the MESA scheme is extended here to 
the nonlinear Euler equations so that the design accuracy required for an efficient direct 
numerical acoustical simulation can be quantified. 

In addition, as the accuracy requirements increase so does the complexity of the non- 
linear simulator. Many subtle errors can be introduced, particularly for nonlinear systems 
which have few analytical solutions available for testing with. Therefore, the Method of 
Manufactured Solutions 13 is applied here and found to be a useful tool for developing and 
testing very high order methods on complex nonlinear systems. 


2. Problem Description 

The following form of the two-dimensional nonlinear Euler equations are used 

( 2 . 1 ) 
(2.2) 

(2.3) 

(2.4) 

This form does not explicitly use the energy primitive variable since the equation of state 
was used to eliminate it, resulting in only four unknown variables: p, p, u, u, which are the 
density, static pressure, x-velocity, and y-velocity respectively. The ratio of specific heats is 
defined as 7 = 7^ = 1.4. 

Ordinarily, the source terms (Qi, Q21 Qs, Qa) are zero. But, for testing purposes, these 
terms will be modified so that the following analytical solution is defined: 


p(xj y, i) = a\ cos(k x nx) cos(A;^7ry) cos(kf 7ri) + c\ (2.5) 

n(x, y , t) = cos(k x nx) cos(k y ny) cos(k t nt) + C2 (2-6) 

v(x , y, t) = <7,3 cos(A; x 7nr) cos(k y ny) cos(kt^t) + C3 (2.7) 

p(x , y, t) = cos(k x 7 Tx) cos(kyTry) cos(k t 7 rt) + C4 (2.8) 


This simple harmonic solution is used to test the resolution of each method as their 
design accuracy is modified and it enables comparisons to be made with previous linearized 
results which also used this test 12 . The source terms, Qi , are then found using computer 
algebra and the resulting code is verified using the Method of Manufactured Solutions 16,17 . 
The source terms corresponding to this analytical solution are shown in Appendix A. 

It is important to insure the density and pressure variables, p(;r,y, t) and p(x,y,t) y are 
greater than zero when choosing the free parameters, (a* and c*), because the numerical 
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approach will otherwise be numerically unstable due to excessive roundoff error from division 
and multiplication by small numbers. The linearized Euler equations were previously solved 
with a linear MESA scheme 12 which did not have this source of instability. However, since 
p and p are not perturbation quantities in Eqs. (2. 1-2.4), they will never be negative in 
practical applications and therefore this instability is only of concern while imposing the 
artificial solution, Eq. (2.5). 

3. Numerical Approach 

The numerical approach described here will produce single step, explicit, O = N(s + 1) — 1 
(V integer s > 0 and N > 1) order of accuracy in space and time methods for the solution 
of the two-dimensional nonlinear Euler equations using a N x N stencil containing the 
primitive variables and their mixed spatial derivatives, a r o...^ y o ...3 ♦ 

The primitive variables and their spatial derivatives are advanced in time on a staggered 
grid 12 such that time advancement is always performed in the center of each stencil by: 


where 


d" +b p(x,y,At) _ y, 1 p , 

01 y k=o K - 

(3.9) 

a“ + V(x,y,A#)u(o:,y,A^) _ ^ 1 rHm 
dx a y b k\ a ' b ’ k 

(3.10) 

d il+b p(x J y,At)v{x,y,At) ^ 1 ^ pv Ajl 

dx a y b h k '- 

(3.11) 

d a+b p{x,y,At) ^ 1 nV Ajfc 

dx“y b £r , 0 kr a ' b ’ k 

(3.12) 

ga+b+ k f{xyt = 0 ) 
a,b,k dx a y b t k 

(3.13) 


The mixed space-time {C[ bk ) derivatives are found by differentiating the governing 
equations in space and time using the Cauchy-Kowalewski technique 18 . The difficulty in 
doing this with nonlinear systems is in differentiating products and this was cited in the 
past as a reason for not using it 19 . Fortunately, it can in fact be performed systematically 
using a multidimensional extension of the Leibniz rule 20 : 


d a + M f{x,y,t)g{x,y,t) 

dx a yH k 


a 





k \ ^-0+(ft-J)+(fc-«=)/(j, y, t) y, t) 

c J dx a ~ l dy b ~i dt k ~ c dx'dyidt 0 


By repeatedly applying the Leibniz rule to the governing equations, the explicit forms of 
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the coefficients, C[ bk , are found. These coefficients are evaluated with the following loop: 
DO k = 1,0; DO b = 0,0; DO a = 0,0 
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DO b : (a + b = 0,1,..., 20) and (0 < a, & < O) 
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END DO 


(3.16) 


After each loop of the k index in Eq. (3.15), the following momentum coefficients are 
known: 

C^ b k with (a, b ) : (a, b = 0, 1, 2, . . . , O) 

C aM witfl (°> 6) : (a, b = 0, 1, 2, . . . , O) 


(3.17) 


But the C% bk and C* bk coefficients are required to evaluate the coefficients, O 0 0 k+1 
at the next k + 1 loop. By repeatedly applying the Leibniz’ rule to Eq. (3.17) an explicit 
expression for the C“ bk and C v a bk coefficients is found, Eq. (3.16), in terms of already 
known variables. Note, however, that these coefficients must be evaluated in the sequence 
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of increasing derivative order (i.e., Find and C^ l)k . in the order (a,b):(a+b=0,l,. . . , 2 

O). By applying this conversion step once for each k-loop increment in Eq. (3.15), we can 
efficiently continue increasing the index k until all the coefficients, Eq (3.13), are found. 


4. Determining the Spatial Derivatives 

The procedure described in the last section will advance the nonlinear Euler equations in 
time at (x,y) with arbitrary accuracy if all the spatial derivatives: 

CL 0 V(«,6) :a,& = O l l 1 ...,0 

are known at (x } y). However, only the spatial derivatives: 

C l M V(«,6) :a,6 = 0,l,...,s 

are available on each grid point, (x uyj). Therefore, the higher order spatial derivatives, 
(a , b > s + 1), need to be interpolated. In addition, since time advancement is performed 
using a staggered grid, both the known derivatives on the grid and the higher order unknown 
derivatives need to be interpolated to the stencil center, (x = 0, y = 0). A simple method 
for accomplishing this on a N — 2 size stencil was shown previously 12 . It is completely 
generalized here to include any size stencil (N > I) using any number of spatial derivatives, 
s > 0, and at any location such as on solid boundaries which require one-sided stencils. 

A one-dimensional osculating spatial interpolant using the primitive variables and their 
spatial derivatives at the grid points, (xq, X \, . . . , :r/v_i), is: 

A r — 1 fc($+l)+s 

_ V" F 

°o,o,o ~ 2s 2s 

k—0 i~k(s~ f-1) 


' k— 1 


5 + 1 


Qu ( (ar-x*r fc(i+1) 

J= 0 / 


(4.18) 


where the are computed as follows: 

DO k = 0, N — 1; DO i — k(s 1), k{s + 1) ■+ sj 
DO j = 0,i - k(s + 1) 

Qu = (i 

END DO] END DO ; END DO 


(4.19) 


DO m = 1,N — 1 ; DO k = 1 , AT — m 

DO i = (s + l)(m + k — 1), (k + m)s + m + k — 1 

DO j — i + 1 — (s + 1 )k, i -f- s + 1 — k(s + 1) (4.20) 


END DO; END DO] END DO] END DO 


The real value of this algebraic form is the simplicity of interpolating the derivatives; 
Ordinarily, by the product rule of differentiation, each differentiation of Eq (4.18), would 
result in an equation that grew in size by a factor of N, up to N° terms, and would render 
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very high order methods intractable. Instead, the following explicit algebraic form avoids 
this issue: 

A 7 — 1 A“(s+1)+5 

Cf'oV = E E QijZ(t,i,k (4-21) 

k= 0 i-fc(s+ 1) 

where the function is independent of space and time and can be computed simply as 

follows: 


DO a. = 0, (N- l){s + l) + s 
DO i = 0, s 

a — 1 

Zajfl = (x - Io) (f_a) n (* - e ) 

c— 0 

END DO 
END DO 

DO a = 0,{N — l)(s + 1) + s 
DO k = 1, N - 1 
DO i = k(s + 1), k(s + 1) + s 
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END DO 
DA^D DO 


1 


2-^p — 0 






This technique is an extension of Hermitian divided-difference interpolation 21 and it 
was derived by repeatedly applying the generalized Leibniz’ rule in Eq. (3.14). For a fixed 
stencil size on a uniform grid, Eq. (4.22) needs to be computed only once since it is constant 
for all space and time. This results in significant computational savings since only the 
variables require further evaluation as the stencil location changes, but even these stay 
constant as higher order derivatives are evaluated using Eq. (4.21) on the same stencil. 

In addition, higher order interpolants may reuse the information from lower order inter- 
polants to efficiently and adaptively increase local accuracy: 


r>fn{x) _ rifN- iW 

^0,0,0 “ u a,0,0 


(AT-l)(s+l)+s 

E 


Qi,iZa,i,k 


i=(AT- l)(s+l) 


(4.23) 


Multidimensional spatial interpolations are efficiently accomplished using the tensor 
product of a series of one-dimensional interpolants 11,22 . 
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5. Explicit forms for the Analytic Solution and its Source Terms 

Computer algebra can calculate the derivatives of the analytical solution, Eq. (2.5), and 
the source terms, C^ 1 bk _ i in Eq. (3.15), and then automatically write the necessary FOR- 
TRAN code 23 ’ 13 , However, for very high order algorithms the code becomes too large and 
complicated for most compilers and therefore it is necessary to develop an explicit form of 
those expressions. 

First define the Lagrangian operator as: 

(n”Jo Mod(a, 4 ) - i) (rTU+1 Mod(a, 4 ) -i) 

L a ri = - 7 ( (5.24) 

(nnu 1 ( n - o) (Tiz-n+i( n - o) 

And then define: 

CT [ * Z) = (Ld,0 COS (ij>z) - L d ,i sin - L d ,i cos {t/>z) + L d , 3 sin(V>z))^ d (5.25) 
C^* z) = (Ld,o sin (ipz) + L d , 1 cos (t()z) - L d ,2 sin (ipz) - L d ,s cos(ipz))ip d (5.26) 


where is any constant, d = a, 6, or, fc, and z = k x irx , k y 7vy, or, ktnt. 

The derivatives of the assumed analytical solution, Eq. (2.5), may then be written as: 


CaTf = taoSboSkoCi + aAk^)^k^)Hk^) k ^ (5.27) 

SaoSboS ko c 2 + a 2 (k x Tr) a (k y ir) b (ktir) k C c a ° (5 28 ) 

3 (A x 7r) a (A:j / 7r) 6 (fc f 7r) fe C“ s(fe7r:E) C , “ s( *»**>£*“(****) (5.29) 

c p a Tk = SaoStnfoc* + ( 5 - 3 °) 


anal 


C'aM = ^aO^boSkoC:) + (13 


where <$ z o is the Kronecker delta. 

The source terms involve products of transcendental functions which may be differenti- 
ated explicitly using the Leibniz’ rule, Eq. (3.14), and are given by: 


n trigi (ipi z)trig 2 (faz) 
° d 


n trigi{^iz)trig2{^2z)trig2{rp3z) 
U d 



ritrigi {rpi z ) n trig- 2 {^2 z) 

U d-i 

) n trigi (rpi z) n trig 2 i$ 2 z) ^trigztyzz) 
U d-i U i-i 2 U i2 



where trigi is either the cos or sin function and ^ is some constant. 

The explicit forms of the derivatives of the source terms using these definitions are listed 
in Appendix A. In this way, only a few pages of FORTRAN code are required to compute 
all the necessary very high order derivatives of the analytical solution and source terms. 

The Method of Manufactured solutions permits any analytical solution to be imposed, 
but note that not all solutions and their induced source terms will be explicitly differentiable 
as they are here. Therefore, additional care is required in the selection of analytical solutions 
for testing very high order methods for otherwise the resulting code will become unwieldy. 
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6. Stability 

Deriving analytical stability conditions for the nonlinear Euler equations is generally diffi- 
cult 24 and necessary stability conditions (but generally not sufficient) are typically derived 
by linearizing the conservative governing equations using their constant coefficient quasi- 
1 inear form: 


OU A dU _ dU n 

at + A & + B % =0 


(6.31) 


where A and B are the flux Jacobians which are held constant for a von Neumann stabil- 
ity analysis. However, in this work, the energy equation is not in conservation form and 
therefore this form of analysis is not available. 

Instead, engineering rules of thumb were applied here and tested numerically. Let the 
Courant number be defined by . 


o = 


'£x_ . 

,Aa- Ay) 


At 


(6.32) 


where c x = rna:v(\\u\\ + c), which is the magnitude of the greatest wave velocity in the x- 
direction over the entire grid (and similarly for c y in the y direction). In addition, for most 
single-step explicit numerical formulations in two-dimensions the Courant number must be 
a < 0.5 to maintain numerical stability 26 . 

Since the Cauchy-Kowalewski procedure is essentially a Taylor series method, a simple 
model problem may be used to demonstrate numerical stability limits. For example 2 ', if 
the test equation (for A < 0): 

y = A y (6.33) 


is time advanced using an n th order Taylor series in time, it will be stable if the following 
series (which is generally convergent as n — > oo) is satisfied: 

||1 + At\ + Q) At 2 A 2 + . . . + Qj) At" A" || < 1 (6.34) 

At low accuracy the maximum stable time step will vary, but it soon approaches a fixed 
value as the accuracy increases. In general, however, high-order explicit schemes usually 
have more severe stability restrictions on At 28 and this is evident in Table 3. 

The following time step size was found to provide numerically stable solutions for up to 
15 fft order accuracy on a uniform Cartesian grid: 


At = 


C 7 


max ; 




where the speed of sound is defined by : 


c= . 7 


pQpjb*) _ / 7(«4 +C4) 

p(x, y, t) \j Cl- fli 


(6.35) 


(6.36) 
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Numerical instability did occur when the solution was under- resolved as shown in Table 
6. Here, the 3 rd order method became unstable by the time t = 799.98, but the 5 th , 7 tf \ 
and 9 th order methods were stable to at least 50 wavelengths of propagation. This suggests 
that nonlinear aliasing effects can be delayed significantly by increasing the accuracy of the 
simulation. Similar results are seen in Table 3. 

In practical applications, it will be necessary to filter out the inevitable unresolved high 
frequencies. However, it is encouraging to observe that despite the very low numerical 
dissipation, filters were not required in this work. 


7. Numerical Results 

A comparison of the effects of design accuracy on the efficiency and resolution of the overall 
simulation was performed by using a 2 x 2 stencil to produce 2s + 1 order schemes. These 
schemes were used to propagate a single wavelength the following distances: 

• A single time step as shown in Tables 1-2, 


• Five wavelengths as shown in Tables 3-5, 


• Fifty wavelengths as shown in Table 6. 


The wavenumbers, amplitudes, and displacements defined in Eq. (2.5) were set to ki = 1 , 
a z = 1 and C{ — 2, respectively except in Table 6 which has c* = 10. 

The expected error, y s , of a method with design accuracy O — 2s + 1 can be extrapolated 
using the data in Table 1 since the design accuracy determines the rate of error change as 
the grid is refined. Starting with the known numerical error, e s , in Table 1 with x s — 2, the 
following relationship: 


Vs ( 2 °(^- 2 )) 


(7.37) 


will provide the expected error at x s = [ Log 2 I • Figure 1 plots this equation. 

Note that, due to computer roundoff, errors less than 10“ 15 are only possible if quadruple 
precision floating point numbers are available. 

By solving for x s in Eq. (7.37), it is possible to estimate the number of grid points 
required to meet a particular error level y s . 


x. = lo ft ~ 1052 Hi + 2 (7.38) 

The total number of grid points in a two-dimensional unit domain required for propagating 
a single wavelength is then given by ( 2 Xs + l) 2 . In addition, each grid point contains the 
four primitive variables p, p, u, and v as well as their, O 2 , spatial derivatives resulting in a 
total memory cost of: 


cost memoT y — 4(2 s + 1) (2 5 + 1) 


(7.39) 
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Table 1. Numerical error at t=0.05, (7=0.25, 4 grid points per wavelength, 25 grid points 
d = 2 = 1 ,ki = 1, At = .05, Single Step 


s 

Order 

Error, e s 

Seconds/Grid Point, c s 

1 

3 

1.5710^ 

.0012 

2 

5 

2.74 10 -4 

.015 

3 

7 

2.54 1(T 6 

.1172 

4 

9 

1.44 10' 8 

.678 

5 

11 

5.40 10~ n 

3 

G 

13 

1.42 10“ 13 

11 

7 

15 

1.78 10“ 15 

36 

8 

17 

3.55 10“ 15 

95 

9 

19 

roundoff 

238 


Table 2. Numerical error at t=0.05, <r=0.5, 8 grid points per wavelength. 81 grid points 
c.i = 2,«j = 1, ht = 1, At = .05, Single Step 


s 

Order 

Error, c s 

Seconds/Grid Point, c s 

1 

3 

9.55 10“ 4 

.0016 

2 

5 

3.9510“® 

.0183 

3 

7 

3.18 10 -8 

.1449 

4 

9 

3.4010“ !1 

.839 

5 

11 

4.85 10“ 14 

3.7 

6 

13 

roundoff 

13 


Table 3. Numerical error at t=10, cr=0.5, 4 grid points per wavelength, 25 grid points 
d = 2, a 2 = 1,/e, = 1, At = 2.4755 10 -2 , 403 Steps, 5 wavelengths 


s 

Order 

Error, e s 

Total Seconds 

Seconds / Grid Point, c s 

1 

3 

UNSTABLE 

UNSTABLE 

UNSTABLE 

2 

5 

1.13 10" 1 

188.22 

7 

3 

7 

1.67 10" 3 

1497.53 

59 

4 

9 

1.50 10“ 5 

8610.51 

344 

5 

11 

5.21 10“ 6 

39125.30 

1565 

6 

13 

UNSTABLE 

UNSTABLE 

UNSTABLE 


Table 4. Numerical error at t=10, <7=0.5, 8 grid points per wavelength, 81 grid points 
a = 2,ch = 1, ki = 1, At = 1.2377 10 -2 , 807 Steps, 5 wavelengths 


s 

Order 

Error, e s 

Total Seconds 

Seconds / Grid Point, c. s 

1 

3 

1.4710“° 

166.90 

2 

2 

5 

3.95 10“ 3 

1392.09 

17 

3 

7 

9.7310“® 

13328.04 

164 

4 

9 

1.6410“* 

64907.56 

801 

5 

11 

1.45 10“ n 

275939.53 

3406 

6 

13 

8.29 10“ 13 

999763.10 

12342 
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Table 5. Numerical error at t=10, <7=0.5, 16 grid points per wavelength, 289 grid points 
a = 2 , «.£ = 1, = 1, A i — G.188 10~ 3 , 1015 Steps, 5 wavelengths 


5 

Order 

Error, e s 

Total Seconds 

Seconds / Grid Point, c s 

1 

3 

1.48 10" 1 

935 

3.2 

2 

5 

2.54 lO -4 

10198 

35 

3 

n 

l 

2.02 10 -7 

80299 

277 

4 

9 

1.29 10~ 10 

4G3537 

1G03 


Table 6. Numerical error at t=100, a— 0.5, 8 grid points per wavelength, 81 grid points 
Ci = 10,aj = l, hi = 1, A t = 5.077 10 -3 , 19C92 Steps, 50 wavelengths 


s 

Order 

Error, e s 

Total Seconds 

Seconds / Grid Point, c 3 

1 

3 

UNSTABLE 

UNSTABLE 

UNSTABLE 

2 

5 

1.05 10“ 3 

G5407 

807 

3 

7 

8.65 10 -6 

40G932 

5023 

4 

9 

1.08 10 _s 

1602013 

20518 


In addition, the total time per grid point for each method, c s is shown in Table 1. However, 
this time reflects the cost of a single time step. As x s increases (grid density) the maximum 
stable time step will decrease for explicit methods (to maintain CFL condition). This is 
accounted for by including the factor, 2 Xs “ 2 , for a total time cost of: 

costume = c s 2 Xs ~ 2 (2 Xs + 1) 2 (7.40) 

Fig. 2 compares the relative costs of propagating a wavelength over a range of error toler- 
ances, 10 _1 ° < y s < 10°. 

Generally, a direct acoustical simulation requires an error of less than 10~ 6 to capture 
both the acoustical and mean flow phenomena. By this standard, the correct propagation of 
a wave cannot be realistically achieved with a 3 rd order method because of the exponential 
growth in its cost as shown in Fig. 2. In addition, additional error will accumulate at each 
time step and as Tables 3-5 indicate, a higher order method is required to propagate out to 
a distance of 5 wavelengths. As Table 6 shows, at least 7 th order accuracy is required for 
propagating 50 wavelengths with 8 grid points per wavelength. 

8. Conclusions 

The procedures showm here provide arbitrary accuracy in both space and time for nonlinear 
Euler simulations using only a few pages of FORTRAN code. The usefulness of these 
procedures hinges on the required accuracy of a given simulation. 

If an accuracy of 10 -10 is required at time t = 0.05 to maintain a given error tolerance at 
some future time, then Fig. 2 indicates the 3 rd and 5 th order methods could not accomplish 
this because of exponential cost growth; Here, the optimal design accuracy is II th order. 
In addition, if the error tolerance is decreased to 10 _1 ° because of the desire to run the 
simulation over a larger time period, it is most optimal to use a or higher order 
method. 
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Figure 1: Resolution Studies 


ORDER 

“♦—3rd 

* 5th 

* 7th 
-^-9th 

- 11th 
— ♦ — 1 3th 
-* 15th 


These results show that a fixed order algorithm approaches its maximal utility and then 
quickly becomes useless. The nonlinear implementation of the MESA schemes described 
here can be used to adjust the design accuracy arbitrarily and it offers an opportunity for 
solving the more demanding acoustical simulations in smooth flows. 
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Appendix A 

Sources terms Q u Q 2 , Qiu and Q 4 in Eqs. (2.1, 2.2, 2.3, 2.4) and their explicit derivatives 
are defined below using the definitions in section 5: 


ryQ: 

— 


r^Q'1 

°uAJfc — 


, U ^ riCOs{k x nx) n cos{k y Ky) ^cos(ktnt) 

(l 1 KfT\K X yt Ot/ 

02^1 k x nk xyt cf s{kl K,) Ky) cf l(fc * - 

t 7 ^cos{ktxt) ^cos(kyny) n shi(k x nx) 

Q’1^2 C(t 

aia2k x nk X ytCl° s ^ kt Kt ^ C^ os ^ ky w) 2 c;sin(2k x irx) _ 

^ t . / ^cos(A-/7rO / ^cos(A* i 7rx) /0 sin(/ry7r2/) 

fllC 3 ay 7 TrCxyf L/ a '-'fc 

a 1 0.3 ky'Kk xyt C c j ? s ^ ktnt ' >2 C^ os(Al ^ x)2 Cj in(frj ' 7ry) 

-a 4 /c, 7 r^ f 6 ;r ( ^" x) C s( ^^V;f l(A ''"' ) - 


G 4 C 2 7T k X yf C a 


'Sm(k x 7rx) £jCOs(k y 7ry) f-ycos(kt'xl) 


cr 


,sin(A* x ?rx) ^ f cos(A' , y ?ry) ^cos(Av ?rf) 


°A- 


/^Q 3 

U aAA" “ 


t 7 ^cos(A~f 7 rt ) 2 ri cos(k y ny ) 2 ^cos^x nx) sii\(k x irx) 
(l 2 (l 4 K x nK X ytL, k U b U a 

T a 2 a 4 A: x Trk X yt C c k ° s(kt ^ C 0 cos( * I,ri) si " { *‘ * j) - 

a 4 C3^7rA: I ^Cr (A ^ x) Cf (fc ^ y) C s(fc ' ,r#) - 

a 3 c 4 7Viw^ fe ” ) ^ n ^^ )c r ( * ,,rf) - 

qofl4 7 7 n cos(k t xt ) 2 n cos{k x xx ) 2 ^\T\{ 2 k v Ky ) 2 

2 0 (J; 

A 03041 • / ^cos(A*f7r/) 2 rt cos(A" x 7rx) 2 rr sin(2A: y 7ri /) 2 

/ y~ 2 K y 7 ^ Kx y fL 'k ° a W 

j 1 ^cos(A* x 7 rx) J ^cos(A , ?/ 7 r 2 /l ^smtktnt) 

-a 2 cik t 7Tk X ytCa K J C b y - 

j , /0 cos(A“ x 7 rx) / - Y cos(A , y 7 ry) / ^sin(A-( 7 rf) 

aiC 2 kt 7 Tk xyt Cu c b c k 

0 , , / ^cos(fc x 7rx) 2 / -YCOs(A* y 7r?/) 2 ^cos(A*f7rq sm(k t nt) 

Z(l\ d 2 Kf 7T K xy f Oa Oju 

j , ^cosflfe* ttO ^cos( A~ y Try) n sm(k x nx) 

d^ruxTr K xy f\y ^ L/ b Oa 

2a 2 c 1 c 2 fc I 7rA: xyt Cr (fc,rt) Cr {fcy ^ ) ^ in(fcl,ri) - 

„ sycos(ktni) siCOs(k y ny) s-ysm(k x nx) 

Qlt^Kx'tK'xyt'^k 

0 2 » 1 / ^cos(^ 7 rf) 2 / oCos(A‘ y 7 r 2 /) 2 ^cos(A , J ; 7 rx)sin(A , x 7 rx) 

K x 7 TK X yt O^. O a 

4«, a 2 c 2 fc x 7r A: x , f C a cos( ** * x) sin( *‘ ^ - 

q 2, 7 ^cos(/p/7r<) 3 ^005(^3; 7T2/) 3 ^- Y cos(A: x 7rx)cos(A‘ r 7rx)sin(A: x 7rx) _ 

oOiO<2^x7Ta^^f O^. O' a 

oscxca^jr^cT^^cr 1 *-^^"^^- 

a 2 Cl C3^7r^ f C, cos(fcfrt) CS os( ^" l) C 6 sin(fc, ' 7r?/) - 
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t 7 s^CQsiktnt) ^cos (fc x ttx) y-ySin(A’ y ttt/) 

aiC2C 3 /^7rA: x ^G* 'C^ G 6 - 

0 , , ^cos(ktTTt ) 2 —cosikxirx ) 2 ^cosikyxy) s\n(k v ny) 

2.(i ( 2G'§C\ fcyff fcxyt^-' fa Ay</ ^5 

O 7 7 ^*cos(Av Trf) 2 ^,cos(A- x 7i\r) 2 y-ycos^y Try) sm(k y 7ry) 

2aia3C2k y nk xy tC k C a G 6 - 

o 1 u ricosiktnt} 2 n cos(k x 7TT) 2 ^cosiky jry) sm(k v 7ry) 

Z(l] a 2 C[] k'y ft KxyfLsfc Is a Ay ^ 


3(7 1 Q^l^kyTT kxyt C fc 


s-ico^kt 7rf ) 3 >-,cos(A , x ttx) 3 y^fCOs(k v Try) co s{h v Try) sin(fc v 7n/) 
Ay/, Ay a Ay i. 


c: 


Q4 _ 


aM,k 


-a$c\i 


j i* y-y€Os(A• x 7 ^Jr) / -tCos(A , I/ 7 ry)^sin(A' , / 7 rf) 
^TT/Cj.y^Aya Ay ^ 


, , ^cosfA^ ttx) fyCos(ky ny) y-ysi n (kt wt) 

Q‘\C'^Kf'ITKxytks(i Ay ^ Ay^ 

o / i ,ocos(fc x nx) 2 f^co^ky ny) 2 n cos(kt nt) sm(k t nt) 
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a 2 c { c z k x nk xyt C™ {k ^ 

1 j fyCos(kt7rt) / ^cos(A , l/ 7ry)^- Y sin(A* x ttx) 

AA 1 C 2 C 3 kx ft G ^ G^ G u — 

0 7 t ^^(A-/ 7 rf) 2 / -tCos(A‘ i/ 7 r^} 2 r ,cos(A , x 7 rx)sm(/r x 7 rj) 

za2 Ci K x ft /C XJ/ ; O ^ Ay t , 

2 (i 1 Q 3 <*2 kx ftk xy tC p, 

2q> 1 U ( i C 3 k x ft k X yf C 


O ^ 7. y-tCOs(A*77Tf) 2 / ^COS(A^7ri/) ^ 

2aiasC2k x ftk X y t C k C b G a 


cr 

2 -cos(fc x 7rx) sm^Tr*) 


rCOs(A*/ 7r0 2 y^cos(A" v ^t/) 2 ^cos(A* x 7rx) sin(A* r ttx) 
Ay;. Ay a 


o 7 7 / oCOs(Ar77r^) 3 / - t cos(A , 3/ 7r7/) 3 / - t cos(A , x 7rx) cos(A' x 7rx)sin(A* x 7rx) 
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